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ABSTRACT 

We extend the existing analytical model of reionization bv lFurlanetto et all (120041) to include the biasing of 
reionization sources and additional absorption by Lyman Limit systems. Our model is, by construction, con- 
sistent with the observed evolution of the galaxy luminosity function at z < 8 and with the observed evolution 
of Ly-a forest at z < 6. We also find that, for a wide range of values for the relative escape fraction that we 
consider reasonable, and which are consistent with the observational constraints on the relative escape fraction 
from lower redshifts, our reionization model is consistent with the WMAP constraint on the Thompson opti- 
cal depth and with the SPT and EDGES constraints on the duration of reionization. We, therefore, conclude 
that it is possible to develop physically realistic models of reionization that are consistent with all existing 
observational constraints. 

Subject headings: cosmology: reionization: theory - methods: analytical 
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1. INTRODUCTION 

Cosmic "Dark Ages" began when the Universe was only 
380,000 years old, at the epoch of recombination, when cos- 
mic plasma and radiation have cooled enough to allow elec- 
trons to combine with protons to form stable neutral hydrogen 
atoms. Lasting for some 200 million years, the Dark Ages en- 
compass time when the dark matter was collapsing into bound 
objects (halos) due to gravity, but no stars were formed yet. 
The epoch of reionization (EoR) started with the formation 
of the first star-forming galaxies and quasars. According to 
current understanding of reionization, the ultra-violet radia- 
tion from these objects was the primary source of energy that 
stripped electrons off the hydrogen atoms, leaving most of in- 
tergalactic matter highly ionized. 

There are three primary observational probes that provide 
information about the EoR. The first one is the Lyman-alpha 
forest, which can trace the distribution of neutral hydrogen 
in the universe and the small-scale structure of ionized bub- 
bles, but only up to z ^ 6 (i.e. only late stages of the EoR can 
be probed with this method). Another probe is the surveys 
of Lyman-alpha emitters, which can statistically constraint 
the size distribution of ionized bubbles from the clustering 
strength of Lyman-alpha-bright galaxies. Finally, Thomson 
scattering of CMB photons off the free electrons produced 
during the reionization, measured by WMAP and SPT, can be 
used to place global, integral constraints on the overall timing 
and duration of the reionization epoch. 

The most promising type of observation, that should be- 
come possible in the near future, is detection of the 21 cm 
transition of neutral hydrogen. The 21 cm emission is a direct 
probe of the large-scale distribution of neutral hydrogen dur- 
ing the reionization era, and is, therefore, a complementary 
probe to the Lyman-alpha forest at lower redshifts. Hence, 
theoretical studies of EoR are timely and relevant. 

The epoch of reionization can be studied theoretically with 
numerical, semi-numerical, and analytical models. Until re- 



Department of Astronomy & Astrophysics, The University of 
Chicago, Chicago, IL 60637 USA; kaurov@uchicago.edu 

^ Particle Astrophysics Center, Fermi National Accelerator Laboratory, 
Batavia, IL 60510, USA; gnedin@fnal.gov 

' Kavli Institute for Cosmological Physics and Enrico Fermi Institute, 
The University of Chicago, Chicago, IL 60637 USA 



cently, numerical simulations of the EoR have been limited to 
small volumes or low resolution (Trac & Gnedin 2011). Be- 
sides that, simulations are way too computationally expensive 
to be used in Markov Chain Monte Carlo (MCMC) analysis 
for obtaining constraints on cosmological parameters or for 
producing multiple mock data sets for the future observations. 
Here semi-numerical and analytical models come into play. 

As the mass in the collapsed objects in the universe was 
growing, the total number of ionizing photons was increas- 
ing too. At some point before z ^ 6, the number of pho- 
tons became sufficient to ionize all of intergalactic hydro- 
gen. Direct calculation of the number of ionizing photons 
and intergalactic baryons allows one to track the global ion- 
ized fraction (c.f. iKuhlen & Fau cher-Giguere 2012, and ref- 
erenced therein). However, all the information about mor- 
phology of ionized regions is completely lost in such straight- 
forward calculations. A more detailed treatment of this pro- 
cess can be achieved with the analytical model based on the 
Press-Schechter like formalism, pioneered by Furlanett o et alJ 
(2004, hereafter FZH04). This model allows one to track the 
bubble size distribution and the mean free path of ionizing 
photons, as well as the power spectrum of ionized regions, 
which can be used for calculating secondary anisotropics of 
the CMB and for the expected 21 cm signal. In this paper we 
present a logical extension of this analytical model that ac- 
counts for the biasing of the ionizing sources and additional 
absorptions by Lyman Limit systems. 

The paper is organized as follows. In ^the brief outline 
of the FZH04 model is presented. In ^we describe the key 
steps of our method, with the details of our implementation 
to follow. In ^the results are presented and we discuss the 
benefits and possible applications of our model. And in ^we 
conclude. 



2. THE FZH04 MODEL 

The FZH04 model is based on the idea of calculating the 
total number of emitted photons in some large enough region 
of space with known over-density and comparing this number 
with the number of baryons in that region. If mcou is the mass 
of the collapsed objects and m is the mass of some region of 
space, the condition for the region to be ionized is formulated 
as: 
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(1) 



We can interpret ( as ratio of the total number of hydrogen 
ionizing photons produced by one collapsed hydrogen atom 
N^/c to the number of photons required to ionize one hydro- 
gen atom Nifii (which would be 1 if there are no recombina- 
tions), 
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In the presence of recombinations. 
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where free is the average recombination time that depends, 
among other factors, on the clumping factor of the gas. Both, 
A^^/e ™d ^(/H can, in principle, vary in space. Notice that the 
gas clumping is only important in the regime when recombi- 
nations dominate over the first ionization. In this paper we 
follow FZH04 and neglect recombinations. 

Now let's consider some region in space, and neglect re- 
combinations for now. It is ionized if the number of hydro- 
gen ionizing photons in it is greater than the number of hy- 
drogen atoms times N^/c- In order to compare these num- 
bers, we first have to find Weou, which is the function of lo- 
cal over-density. It can be done using, for example, Press- 
Schechter m odel (Pr ess & SchechteiilI974tlBond et alj|1991l: 
lLacev& Cole.. 1993.). 



/c 



coll 



: erfc 



(2) 



where a^{m) is the variance of density fluctuations on the 
scale m and M^i„ is the minimum mass of an ionizing source 
(we discuss this mass later in ^3.2l FI. Hereafter we chose con- 
vention (adopted by FZH04) in which a(m) is evaluated at 
the reference redshift z* = and Sdz) is a function of redshift. 
Such convention is convenient because Sdz) becomes the only 
quantity that explicitly depends on time. 

Hence, the condition that the bubble of mass m at redshift z 
with over-density 6m is ionized becomes 



^ < /coll- 
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Equation (|3]l can be rewritten as a constraint on local den- 
sity: 

S,„ > Sdm,z) = Sdz)-V2K(0 [aHM^i,)-a\m)] ''^ , (4) 

where K.(Q = erF' (1 - <^~'). Using the excursion set formal- 
ism and the function 5.c(m,z) as a barrier, one can find the 
distribution of the masses (and scales) of ionized regions. 

More detailed description of this method can be found in 
FZH04. In the rest of this paper, we assume that the reader is 

^ Hereafter we adopt the FZH04 notation of using a lower-case m to denote 
the mass of an ionized region; we also use a standard convention of adopting 
an upper-case M for the masses of dark matter halos and galaxies. That makes 
our notation somewhat non-intuitive, but we consider that preferable to using 
a more logical but totally unfamiliar notation of labeling the masses of ionized 
regions with M and masses of galaxies with m. 



familiar with the excursion set formalism, the concept of bar- 
rier crossing, and Monte-Carlo sampling of possible random- 
walk trajectories as a method to compute the first barrier 



crossmg. 



3. OUR APPROACH 



The method described in the original work by 
iFurlanetto et al.l (120041) assumes the constant photon per 
collapsed baryon ratio. Later in Furlanetto et al. (2006) it 
was extended by making ( a function of halo mass. However, 
the model based on excursion set formalism allows one to 
take into account the efficiency of photon production as a 
function of three parameters: redshift, local over-density, and 
the bubble size. Below we give a general overview of our 
approach step by step; it is followed by detailed explanation 
of each component. 

1. In the first step we find the UV luminosity-to-mass ra- 
tio of galaxies (i.e. ionizing sources) as a function of redshift. 
One can do that by using abundance matching between the 
halo mass function (Tinker et al. 2008) and the observ ed lu- 
mino s ity fun ctions of high redshift galaxies ( Bouwens et al.l 
120071 12010I) . The main uncertainty in this step is the ex- 
trapolation of the abundance-matched luminosity-to-mass ra- 
tio to higher redshifts, since observations only cover the red- 
shif t ran ge z < 8, while reionization starts at earlier times. 
In ^3.11 we describe in details the abundance matching pro- 
cess and our method of extrapolation to higher redshifts that 
is based only on the theoretically computable evolution of the 
halo mass function, without any additional parameters. 

2. The observations provide us with the information about 
the galaxy UV luminosity at the wavelength ^ 1600A. But 
the ionizing radiation has the wavelengths below 912A. To 
convert one into another we need to know the typical spec- 
trum of the galaxy and the relative escape fraction of the pho- 
tons at these two wavelengths. For the typical spectrum we 
used Starburst99 ( Leitherer et al. 1999; Vazquez & Leithere^ 
I200I iLeitherer etal.l 120101) and found that ratio of fluxes at 
A < 912A and A = 1600A is r,„, = 0.241 (see Appendix for 
details). For the escape fraction, w e use a simple, two param- 
eters model which is discussed in ^3.21 Combining these two 
factors we can rewrite luminosity as: 

^<912a(^.^) = '•int/csc,rel(M)Li6oo^(M,z) 

3. In this step we calculate the number of ionizing pho- 
tons. In FZH04 the authors used the parameter (, which 
measures the ratio of ionizing photons to collapsed hydrogen 
atoms. We find it more convenient to use another parameter, 
N^/H, the ratio of the number of hydrogen ionizing photons 
to the total number of hydrogen atoms in the region (these 
two parametrization are related to each other by the collapsed 
fraction). The luminosity-to-mass ratio for ionizing sources, 
calculated in the first step, allows us to find N^/a at the mean 
cosmic density as a function of redshift. 



N^/h(z)-- 



J L^<)n(M,z)^(z)dM 

"H.O 



(5) 



where dn/dM is the halo mass function and where hh.o is the 
mean number density of hydrogen atoms. 
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At this point we need to account for the halo bias. In over- 
dense as well as in under-dense regions the halo mass func- 
tion is different, which makes N-f/H ^ function of local over- 
density 5. Details of how the bias is implemented in our cal- 
culations are provided in 33.31 Another important factor is 
the absorption of ionizing radiation by Lyman Limit systems, 
which are not modeled in the FZH04 formalism and need to 
be included as a separate component. This absorption de- 
pends on the mean free path of an ionizing photon at a given 
redshift and the size distribution of ionized bubbles. In 33.41 
we describe in details what model for Lyman Limit systems 
we use and how we calculate the fraction fLLsiz,R) of ioniz- 
ing photons that can reach the boundary of an ionized bubble 
of size R (in comoving units) and hit a neutral hydrogen atom. 
Taking the bias and Lyman Limit systems effects into account. 
Equation (|5]l becomes 



dn 



lL^9nA(M,z) — (z,S,R)dM 
nH,o(l + ^) 



(7) 
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where (z, S, R) is a halo mass function at overdensity 6 on 

dM 
scale R. 

4. In the most general case N^f^ is the function of redshift, 
local over-density, and scale (bubble size). The bubble of size 
R at some redshift z with some over-density S is ionized when 
Nj/uiz,S,R) reaches the fraction of non-collapsed hydrogen 
atoms. 



Fig. 1. — Illustration of using the abundance matching to derive the 
luminosity-mass relation. Cumulative mass function of dark matter halos (top 
panel) and the cumulative UV luminosity function of galaxies (right panel) 
can be matched to obtain the relationship between the UV absolute magnitude 
and the halo mass (central panel). 

sic (i.e. unabsorbed) spectrum of the galaxy L'J;" is known. 



L{< 912A) 



N^,n(z,S,R)=l-f,, 



= /e.«912A)^^!:^ 



where /cou is taken from Equation (|2| with Mjj,i„ being the 
minimum mass of a halo that is capable of r etaining its gas 
after ( local) reionization jHoeft et al.l 120061; lOkamoto et al.l 
I2008al) . In this paper we take that mass to be IO^Mq, since 
the minimum mass of a halo capable of retaining its gas has 
not yet been calibrated as a function of local conditions; how- 
ever, in the future the model can be improved by adopting it 
as function of redshift and other parameters. This correction 
is not particularly important, though, since /coii is expected to 
be small. 

5. The last step in our model is computing N^/h- In the 
original FZH04 model ( was constant and N^/n was a func- 
tion of z only , hence Equation (Q could be integrated directly. 
In the case when Nj/u depends on R (or, equivalently, on m), 
6, and z, it is no longer valid to directly integrate Equation ^ 
because both 6 and R/m can be functions of time. In addition, 
ionized bubbles can merge, thus combining two separate bub- 
bles with nil and m2 into a single bubble of mass m = mi+mz- 

In order to account for both of these effects, we need to find 
the barrier at time f + Af from the information about the bubble 
distribution at time t. We describe our approach in ' 



3.1. Luminosity-Mass Relation 

For a given galaxy, the production rate of ionizing photons 
can be derived from the observed UV luminosity if the intrin- 



= /e.«912A)-^"^<'^^2^>^^'"'(^^«0^> 



= /esc« 912A) 



AL-'(1600A) 
o ^ L''"( < 9 1 2 A) ALf ^ ( 1 600 A) 



AL'"'(1600A) /esc(1600A) (£) 



= /esc,rel-7^ALf (1600A), 



(8) 



where is the average energy of a hydrogen ion- 

izing photon, Tint = {U'^^^^.J \Uf{\6Q0k) is the intrin- 
sic ratio of the ionizing to UV luminosity, and /esc.iei = 
/esc(< 9 12A) //esc (1 600 A) is the relative escape fraction (c.f. 
iGnedinet al]|2008al) . 

We discuss our model for the relative escape fraction be- 
low, in 33.21 For our choice of the intrinsic galaxy spectrum 
(described in the Appendix) r-mi = 0.241 and (£')io„ = L93Ry. 

At redshifts where the galaxy luminosity function is actu- 
ally measured, all we need to do is to integrate Equation (O 
over all galaxies. However, in Equation ^ we need to know 
My at all redshifts, between the observational windows and 
beyond the highest redshift (z « 8) at which the reliable ob- 
servational determination of the galaxy luminosity function 
exists. 

In order to extrapolate beyond z w 8, one can use several 
approaches. One can, for example, fit the observed evolution 
of the galaxy luminosity function with some parametric form 
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and extrapolate that parametric fo rm mathematically, like in 
iKuhlen & Faucher-Giguer j ( l2012h . A limitation of such ap- 
proach is that any extrapolation depends on the exact mathe- 
matical form of the fitting function; different fitting functions 
will produce different extrapolations, thus biasing the result 
in an unpredictable way. 

It would be preferable to base the extrapolation on a physi- 
cal model. In the Press-Schechter framework, the only phys- 
ical component is the evolution of the halo mass function. 
Thus, a physically-based method of extrapolation is to assign 
galaxies to dark matter halos and use the evolving halo mass 
function to compute the evolving galaxy luminosity function. 

To find the relation between the halo mass and the galaxy 
UV luminosity one can adopt a widely used method of abun- 
dance matching: given the cumulative halo mass function and 
the cumulative galaxy luminosity functions, halos of a par- 
ticular mass M (for some mass definition, not necessarily the 
virial mass) can be identified with galaxies of some luminos- 
ity that have the same number density. This matching does 
not have to be fully deterministic - scatter can be included in 
the relation as long as we treat the abundance-matched lumi- 
nosity L{M) as the average luminosity of all galaxies that live 
in halos of mass M. 

Such an approach has been used, fo r example, by ' GnedinI 
( l2008i) and lVolonteri & GnedinI ( 120091) . and we show an up- 
dated version of abundance matching betwee n the galaxy UV 
luminosity and the halo virial mass in Figure lTTI for four val- 
ues of cosmological redshift. 

However, as is apparent from Figure 13.11 the luminosity- 
mass relations obtained by abundance matching are redshift- 
dependent. Hence, any extension to higher redshifts will re- 
quire extrapolation and will again be subject to unknown bi- 
ases; 

iTrenti et aP (|2010) circumvented that limitation by match- 
ing luminosity functions not with the actual halo mass func- 
tions, but with the difference between the two halo mass func- 
tions at two redshifts separated by a fixed time interval of 
200 Myr. With such matching, the luminosity-mass relations 
become redshift-independent. 

In this paper we propose another method of extrapolation. 
Our main motivation is to find a physically plausible relation 
between the galaxy luminosity and some property of a dark 
matter halo, and then match that property to the halo mass. 
One such property may be mass within a fixed proper density 
threshold. The cosmic overdensity A changes simply due to 
the expansion of the universe, while stars presumably form 
in the interstellar medium of galaxies at densities that do not 
depend on cosmic expansion. 

Hence, if Ma is the mass within the radius enclosing the 
cosmic overdensity A (so that M200 is the virial mass at high 
redshift), we can obtain a relation between the galaxy lumi- 
nosity L and Ma by abundance matching, 

n{> L) = n{> Ma), 
and choose A to correspond to a given physical density Aphys, 

For our fiducial value we choose Aphys = 6.9 x 10^ («;, w 
0.2 cm"^), which corresponds to the cosmic overdensity A = 
3200 at z = 5 and overdensity A = 200 z = 14. 

With that choice of the matched halo mass, the luminosity- 
mass relation becomes redshift independent - i.e., if we use 



^Api,y, in Figure [TT] instead of M200, all four lines fall on top 
of each other (we do not show such a figure since it is trivially 
redundant). 

Given the relation between L and Ma, we can derive 
L(Mvir,z) at any redshift by abundance-matching of two the- 
oretical mass functions, n(> Ma) and n(> Myii), which are 
known at any redshift, i.e. 

L(M™,z) = L(MA(Mvir,z)). 

This procedure does not involve any extrapolation in time, 
and, hence, all time evolution comes through the physical evo- 
lution of the halo mass function. 

3.2. Escape Fraction 

The most sensitive parameter for any reionization model is 
the escape fraction of ionizing photons from galaxies. Dif- 
ferent definitions of f^^c exist in the literature. When the ac- 
tual ionizing emissivity of all stars in a galaxy is known, f^^c 
stands for a fraction of the radiation that leaves the galaxy 
- it is called the absolute escape fraction then. In reality, 
the total ionizing flux from stars in high redshift galaxies is 
unknown. The observed quantity is the UV luminosity (at 
wavelength A ~ 1600 A) of a galaxy, and, hence, we are inter- 
ested not in the absolute, but in the relative escape fraction, 

/escrei = /esc(< 9 12A)//esc( 1600A) that enters Equation 

A large body of work exists on modeling the es- 
cape fractions from high r edshif t galaxies (see, e.g. 
iRazoumov & Sommer-Larsen' '2010*; 'Srbinovsky & Wyithd 
2010; Wyitheetal. 2010; Fernandez & Shull 201 1 
IMitra et al.i 12012 , for the latest models). In this paper 
we adopt the second-to-the-simplest model with only two 
parameters - the minimum mass M^m and the amplitude 

/esc.iel.O- 

/esc,el(M)-|Q^ M<Me,,. 

The mass M^rit in Equation (|9]l may play multiple roles: it 
can be the critical mass Merit ^ 1O"M0 from (iGnedin et alJ 
l2008bh who found that dwarf galaxies below that mass have 
very low escape fractions; alternatively, it can be the minimal 
mass of FZH04 that corresponds to the virial mass of a halo 
with temperature 10"* K (the critical temperature for produc- 
tion of ionizing photons. IChoudhurv et alj|2008l) . or it can be 
a mi nimum mass of a halo capable of retaining photoionized 
gas (lOkamoto et alJl2008bh . 

3.3. Halo bias 

The bias factor modifies the halo mass function in the re- 
gions with non-zero cosmic over-density 5, and, therefore, 
forces the photon production rate to become a function of 
local over-density. To calculate the total luminosity in a 
given region, we need to integrate the luminosity-mass re- 
lation L^gj2^(M, z) from 33.21 over the halo mass function 
to find the average UV luminosity produced by a collapsed 
baryon at redshift z, 

i^<9i2k(M,z)x^(z,6,)dM, (10) 

where we use the bubble mass m instead of its R as a vari- 
able, as a more convenient variable (in this formulation there 
is no need to assume a particular geometric shape for a bub- 
ble). We adopt the mass function from .Tinker et al, (.20081) 



5 



with the halo bias model from [Tinker et alj ( 120101) . modified 
in the following way to ensure the exact conservation of ion- 
izing photons: 



dn 
dM 



dnxinker 



/■Tinker 
/coll 



dM 



(11) 



The factor /^^y//^']"'"^'" is required, since the original Press- 
Schechter form for the halo mass function allows to maintain 
the exact photon conservation independently of the smoothing 
scale by virtue of the property 

(/c'oll(^,))=/c'oll('5 = 0,=0), 



while fits provided in [Tinker et alJ (120101) do not guarantee 
such normalization. We further discuss the conservation of 

photons in the FZH04 model in the Appe ndix. 

The Lagrangian bias that is used in [Tinker et al.l (12010b 
works only on large scales (small ) and small over-densities. 
In ideal case we should know local non-linear Eulerian bias. 
Such information might be extracted from large sc ale simu- 
lations, like it was done in iRoth & Porcianil (1201 Ih . Unfor- 
tunately, at the moment these models can't cover all the red- 
shifts and scales that are considered in our reionization model. 

3.4. Lyman Limit Systems 

In the original FZH04 model the universe remains 100% 
ionized after the end of reionization. The mean free path of 
an ionized photon in such a universe would be limited by the 
cosmic horizon only. In reality, however, the mean free path 
for ionizing radiation is much shorter than t he cosmic horizion 
even at lower redshifts ("Songa ila & Cowi3l201(^') : it is limited 
by absorption in the Lyman Limit systems. Hence, the Lyman 
Limit systems are not modeled by the original FZH04 model 
and need to be accounted for separately. 

If the mean free path is much larger than the size of an ion- 
ized bubble, the majority of photons produced inside a bubble 
will be able to reach the bubble edge and ionize fresh neutral 
hydrogen outside the bubble. In the opposite case, when the 
mean free path is small compared to the bubble size, the prob- 
ability of absorption by Lyman Limit systems is much higher 
and only galaxies near the edge of the bubble contribute to 
ionizing fresh neutral hydrogen beyond the edge. In Appendix 
we calculate the fraction fLLs{z,R) of ionizing photons emit- 
ted inside the bubble that are available for ionizing fresh neu- 
tral hydrogen (and, hence, contributing to N^/u(z,6,R)). For 
such a calculation, distributions of galaxies and Lyman Limit 
systems inside every bubble is needed; such distributions are 
not part of the Press-Schechter formalism and, hence, have 
to be added as external assumptions in the model. In this 
paper we consider two limiting cases: a homogeneous dis- 
tribution (galaxies are distributed uniformly inside a bubble) 
and a highly clumped distribution (all galaxies are located at 
the bubble center). Lyman Limit systems in both cases are 
assumed to be uniformly distributed inside the bubble. 

The average abundance of Lyman Limit systems is not a 
free parameter, h owever; it has been constr ained observation- 
ally up to z « 6 dSongaila & Cowid 120101) . For our purpose 
it is more convenient to use the mean free path due to Lyman 
Limit system, /u s, directl y, hence we adopt Equation (8) from 
ISongaila & Cowig dloTol) : 



/lls(z) = 50 



1+z 
4.5 



-4.44 




(m) 



Mpc. 



Fig. 2. — Illustration of a bairier at time ; (solid line) and two random 
walks (dashed and dotted lines) started at some point (square) which belongs 
to the hairier at time t + At. The dashed random walk crosses the baiTier and, 
therefore, corresponds to the region of space that was already ionized at time 
t. The dotted random walk corresponds to the region which was ionized in 
the time interval between t and f + Af. 

This fit is made in the redshift interval 0-6. In our model 
we extrapolate this fit to higher redshifts, where its accuracy 
cannot be verified, but, at present, there is no alternative. 

To accurately calculate the effective fraction of photons due 
to Lyman Limit systems at some redshift we should know not 
only the bubble sizes but also the merging statistics of these 
bubbles. For instance, if we have a large number of small 
bubbles that merge into larger ones only during the late stages 
of reionization, the effective fraction will be always close to 
unity. In the opposite case, if all galaxies are highly clustered 
at the centers of a few large bubbles, the bubble size will be 
limited by the mean free path and a large fraction of ionizing 
photons will be lost in Lyman Limit systems. Luckily, the 
FZH04 model contains the information about bubble merging. 

3.5. Merging of bubbles 

The procedure of searching the pr ogenitors using excursion 
set formalism was first described in lLacey & Cold (Il993i) for 
the Press-Schechter model. Here we follow the same idea. It 
was al ready applied to the FZH04 model in iFurlanetto & OhI 
(12005b . where the distribution of bubble progenitors was 
found. However, the authors did not use this information for 
correcting the number of ionizing photons. 

Let's imagine that at some time t we have a distribution of 
ionized bubbles in space and assume that there is a barrier that 
can describe this distribution 0. During some time Af each of 
the bubbles was growing according to the photon production 
rate from Equation as well as merging with other bubbles. 
In order to calculate A^^/h for a bubble of size m at f + Af we 
need to know its progenitors, which is illustrated in Figure |2] 
A solid line corresponds to the barrier at time t. If we start a 

' This assumption is ciitical in this theory. 




9.5 10.0 10.5 11.0 

Z 

Fig. 3. — Three reionization histories for a representative model for the es- 
cape fraction (Merit = 3 X IO'^Mq and /csc.rcl.o = 0.4). Solid line - without 
LLS; dashed line - with assumption that galaxies are uniformly distributed 
inside ionized bubbles; dash-dot line - with assuinption that galaxies are 
clumped in the center of ionized regions and LLS are uniformly distributed. 

random walk at some point below the barrier (a black square) 
assuming that this point belongs to the new barrier at t + lS.t, 
the trajectories that cross the barrier at time t will correspond 
to the progenitors of the current bubble. Two random walks 
are shown in the figure: the dashed one crosses the barrier 
and therefore corresponds to the ionized bubble at time f ; the 
dotted random walk does not cross the barrier, thus it corre- 
sponds to the region that was neutral at time f and got ionized 
during the time interval between t and t + iSt. 

If we know the distribution fUm)) of the barrier crossings 
at time t, we can write N^/n as: 

A^.^/h(/ + At,S',m') X m'(l -I- S') = 

Cr^/(('«))xm(l + <5)x 

{N^/iiit,S,m)+N^/nit,5,m)At) dim), (12) 

where we took into account the initial numbers of photons in 
each bubble and their growth during At. The function fUm)) 
in Equation (fT2] i depends on the barrier at time t. Solving 
Equation ( fT2] i for every m' gives us the barrier at time t + At. 

The described approach makes the model more computa- 
tionally expensive than the original FZH04 approach, but still 
much faster than a full numerical simulations. 

4. RESULTS 

The significance of Lyman Limit systems is demonstrated 
in Figure [3] where we show three reionization histories with 
Lyman Limit systems effects turned on and off. The corre- 
sponding barriers at z = 9 are shown in the top panel of Fig- 
ure |4] In the bottom panel the first crossing time distribution 
for the barriers are presented. Additional absorption by Ly- 
man Limit systems leads to the upturn of the barrier on large 
scales (or small (m)), because in large ionized bubbles only 
galaxies near the edge contribute to the further expansion of 
a bubble. Lyman Limit systems affect the reionization history 




a (m) 



Fig. 4. — Barriers (top panel) and first crossing distributions (bottom panel) 
at z = 9 for the three models shown in Fig. [5] 
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Fig. 5. — End of reionization (measured as 99% ionization level, solid 
lines) the duration of reionization (redshift interval between 20% and 99% 
level of ionization, dashed lines) in the parameter space (M^it and /csc.rcl.o) 
of our reionization model. 

at later stages, when the average bubble size is large, signifi- 
cantly slowing down the reionization process and making the 
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f esc, rel, 

Fig. 6. — Observational constraints in the parameter space (M^nt and 
/escrcl.o) of our reionization model (wliicli, by construction, matches the 
observed evolution of the galaxy UV luminosity function for z < 8 and 
the observed opacity of the hy-a forest at z < 6). Solid, dot-dashed, and 
dashed lines show the WMAP+SPT constraint on the Thompson optical 
depth I Keisler et al. 2011) for models without LLS, with LLS and uniformly 
galaxy distribution, and with LSS and highly clustered galaxy distribution re- 
spectively. The shaded region is the SPT constraint on the duration of reion- 
ization from Zahn et al. 1 2012). For reference, we also show with the dot- 
ted horizontal lines the observed constraint on the relative escape fraction at 
z ~ 3 from Nestor et al. 1 2012) an d with the vertical dotted Hnes the constraint 
on the critical mass from Fernandez-Soto et al. (2003) (see text for detail). 
Con straints on the duration of reionization from EDGES ( Bowman & Rogers 
I2012i) are, presently, too weak to be able to exclude any portion of the dis- 
played parameter space. 

end of reionization more gradual. 

The halo bias modifies the halo mass function, resulting in 
higher abundance of massive halos in over-dense regions, and, 
consequently, a higher photon production rate. That lowers 
the barrier on all scales. However, the effect caused by bias 
is minor compared to the one due to Lyman Limit systems. 
Note that the bias effect, in contrast to Lyman Limit systems, 
does not change the total number of ionizing photons, and, 
therefore, does not affect the average history of reionization. 

The computational efficiency of this method allows us to 
explore the parameter space of our model for the escape frac- 
tion (Equation |9l). The results are presented in the Figure |5] 
Solid contour lines show the end of reionization for the model 
without Lyman Limit systems, since the end of reionization 
in models with Lyman Limit systems is less defined (because 
the end of reionization is gradual). Dashed contour lines show 
the corresponding duration of reionization, the redshift inter- 
val between 20 % and 99% levels of ionization (following the 
definition from Zahn et alj ( l2012b ). 

The duration of reionization (Az) varies from 0.4 to L6 
across the parameter space presented in the figures. Current 
constraints on Az from CMB are too broad and can not rule 
out any of our model parameters. The lower limit on Az made 
by EDGES (Bowman & Rogers 2012) and is Az > 0.06 with 
95% confidence level, which is much smaller than any of our 
predictions. 



Other observational constraints in our parameter space 
(Merit and /escrei.o) are shown in Figure|6] For all our reioniza- 
tion models consistent with the WMAP valu e for the Tho mp- 
son optical depth (r = 0.085 ± 0.014 Keisler etin|20lll), the 
measurement of the secondary anisotropies by South Pole 
Telescope (SPT Zahn et al. 2012) does not yet provide any 
significant constraint, because reionization proceeds faster 
than the SPT upper limit on the duration of reionization. We 
also show, mostly for the sake of illustration, in Figure |6] ob- 
servational constraints on out two main parameters (Merit and 
/e,sc.rei,o) obtained at lower redshifts z ^ 3: the actual mea- 
surement o f /escrei.o bvlNestor e t al. (2012) and the constraint 
on Merit by Fernandez-Soto et al., (2003, if we interpret their 
non-detection of the escape of ionizing radiation for the low- 
luminosity sub-sample as a cutoff in the escape fraction for 
low mass galaxies). We reemphasize that the latter two con- 
straints are obtained at z ^ 3 and, therefore, may not apply to 
galaxies at z > 6. 



5. CONCLUSIONS 

In this paper we extend the analytical model for cosmic 
reionization from FZH04 to account for the absorption of ion- 
izing radiation by Lyman Limit systems, for the halo bias ef- 
fect, and for merging of ionized bubbles. 

The halo bias doesn't affect the average reionization his- 
tory since it doesn't change the total number of galaxies (and, 
hence, the total number of ionizing photons). However, the 
halo bias causes a redistribution of ionizing photons between 
over-dense and under-dense regions, and, therefore, modi- 
fies the size distribution of ionized bubbles by increasing the 
abundance of large bubbles and suppressing small ones. 

Additional absorption of ionizing radiation by Lyman Limit 
systems is significant only during the late stages of reioniza- 
tion, when sufficiently large bubbles (with sizes exceeding the 
photon mean free path due to Lyman Limit systems) are abun- 
dant. 

Our model is, by construction, consistent with the observed 
evolution of the galaxy luminosity function at z < 8 and with 
the observed evolution of Ly-a forest at z < 6. We also find 
that, for a wide range of values for the relative escape fraction 
that we consider reasonable, and which are consistent with 
the observational constraints on the relative escape fraction 
from lower redshifts, our reionization models are consistent 
with the WMAP constraint on the Thompson optical depth 
and with the SPT and EDGES constraints on the duration of 
reionization. We, therefore, conclude that it is possible to de- 
velop physically realistic models of reionization that are con- 
sistent with all existing observational constraints. 



This work was done with significant usage of CosmoloPy 
Python packag43- 
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supported by grants from Fermilab, Kavli Institute for Cos- 
mological Physics, and the University of Chicago. This work 
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and arXiv . org preprint server. 
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APPENDIX 
GALAXY SPECTRUM 



In this appendix we compute the fraction photons produced by galaxy that is energetic enough to ion ize hydroge n. We use the 
simpl e fit to the stellar spectrum produced by Starburst99. The Equation 2 in lGnedin & Hollonl(l2012h (Figure 4 in lRicotti et alJ 
l2002h : 



hv 



(Al) 



{5.5, x<\ 

x-'\ l<x<2.5 

0.4x-'^ 2.5<x<4 ^^2) 

2 X 10-V/(exp(x/1.4)- 1), 4 < x 

Since the observations of galaxies (Bouwens et al. 2007, 2010) provide us with information about the luminosity at wavelength 
A = 1600A, we want to convert it into ionization luminosity (wavelengths < 912A): 



L7« 912) = az^Lr(1600), (A3) 

where a is the constant that converts observed value i/L'^**(1600A) into total flux L'"'(< 912A). Let's find a. We know that 
^ s^. So: 



/oo 



from which we find: 

a = 0.241. (A4) 
The average energy of ionizing photon can be calculated as 

{E) = = 1. 93 Ry, ( A5) 

and the emission rate of ionizing photons is 

Nj = aiyL°J"(l600A)/ {E) . (A6) 



* http://roban.github.com/CosmoloPy/ 
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FRACTION OF ESCAPED PHOTONS 

Here we calculate the fraction of photons that end up in Lyman Limit systems and the fraction that successfully reach the 
boundary of an ionized bubble and effectively excite the baryon on the edge. 

Firstly we consider the case when galaxies and LLS are distributed uniformly inside the bubbles. When the mean free path 
(MFP) is larger than bubble size the number of effective photons is proportional to the volume. Another regime is when MFP is 
much smaller, in this case the number of photons is proportional to its surface area. 

If the galaxies and Lyman Limit systems are distributed uniformly in the ionized sphere, what fraction photons can reach the 
surface of the bubble? Assume that MFP at that epoch is /lls- 

Consider a bubble of size R and a galaxy located at distance h from the center. The fraction of photons that wiU reach the 
boundary of the bubble will be: 

^/h^ + R^-2Rhcos(e) = x, 

1 

Uf(h) = j d62TT sm(9) xx^x exp(-x/ /lls)- 
Now we have to find the average over volume: 

/lls(«/A)=^/ AWUn{h)dh. 
In result we have the /lls as a function of /?//lls- Notice, that /lls is the function of redshift. 

The second case is when galaxies are clumped in the centers of bubbles, but LLS are still distributed uniformly. In this case the 
fraction of not-absorbed photons is simple: 

fLLs(R/lLLs) = 1 -eXp(-/;//LLs) 
CONSERVATION OF PHOTONS 

In the original FZH04 model the number of photons is not conserved; i.e. the ionized fraction of the universe calculated with 
the excursion set approach is systemically below the ionized fraction calculated directly, by counting all emitted photons (i.e. 
C/coii(f^ = 0)). This non-conservation is due to trajectories that never cross the barrier - these trajectories correspond to spatial 
locations that are never fully ionized internally, but some ionizing photons are still emitted inside them (there are just not enough 
of them). These photons are not accounted for in the original FZH04 model. To correct for that non-conservation, we extended 
the original FZH04 barrier from Equation (HI with an additional vertical barrier at cr^(M,ni„ + e), where e is small (the barrier is 
undefined if e = 0). Each trajectory that crosses the vertical part of the barrier has some collapsed fraction at this smoothing scale 
and, therefore, the ionized fraction of corresponding region is C/coii(o'^(Mniin + e)), which is less then unity but larger than zero. 
Taking into account all these regions, we restore the exact photon conservation in the original FZH04 model and in our extension 
to it. 



